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ABSTRACT 

Front  tracking  is  an  adaptive  computational  method  in  which  a  lower 
dimensional  moving  grid  is  fitted  to  and  follows  the  dynamical  evolution  of 
distinguished  waves  in  a  fluid  flow.  The  method  takes  advantage  of 
known  analytic  solutions,  derived  from  the  Rankine-Hugoniot  relations, 
for  idealized  discontinuities.  In  this  paper  the  method  is  applied  to  the 
Euler  equations  describing  compressible  gas  dynamics.  The  main  thrust 
here  is  validation  of  the  front  tracking  method:  we  present  results  on  a 
series  of  test  problems  for  which  comparison  answers  can  be  obtained  by 
independent  methods. 
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1.  Introdncdon 

Front  tracking  is  an  adaptive  computational  method  for  modeling  flmd  flow.  Its 
goal,  shared  with  adaptive  methods  in  general,  is  to  obtain  increased  resolution  by  using 
special  computational  degrees  of  freedom  that  (a)  arc  placed  (in  space  and  time)  where 
they  arc  most  needed  and  (b)  fit  the  nature  of  the  solution  as  closely  as  possible.  Within 
this  general  framework,  tracking  is  distinguished  by  the  choice  of  a  lower  dimensional 
adaptive  grid,  called  the  front  or  the  interface,  as  its  special  computational  degree  of  free- 
dom. Thus  for  flow  problems  in  two  space  dimensions,  tracking  employs  a  moving  one- 
dimensional  grid,  i.e.  a  system  of  curves,  in  addition  to  a  two-dimensional  grid.  Further- 
more these  curves  arc  not  purely  geometrical,  but  are  associated  with  physical  waves  in 
the  solution.  They  are  defined  impUdtly  by  the  solution,  and  they  evolve  dynamically 
with  it. 

The  problems  for  which  tracking  is  an  attractive  method  arc  those  containing  discon- 
tinuities and  other  singularities  concentrated  on  surfaces  (curves  in  two  dimensions).  Such 
singularities  abound  in  compressible  fluid  dynamics:  they  include  shock  waves,  contact 
discontinuities,  material  interfaces,  phase  boundaries,  slip  lines,  and  chf-mical  reaction 
fronts. 

What  about  the  internal  structure  of  the  discontinuity?  Sharp  discontinuities  exist  as 
solutions  of  mathematical  equations  (e.g.  the  Euler  equations)  that  describe  idealized  phy- 
sics. The  missing  physics  and  corrected  equations  are  usually  known  (e.g.  the  Navier- 
Stokes  equations),  and  so  are  the  parameters  that  govern  the  correction  (e.g.  the  viscos- 
ity). Such  a  parameter  often  defines  a  length  scale  that  measures  the  width  of  the  discon- 
tinuity. For  many  flow  problems  this  length  is  much  smaller  than  realistic  computational 
mesh  lengths,  so  that  two  possibilities  arise:  either  the  physics  at  subgrid  levels  can  be 
ignored  and  tracking  based  on  the  the  idealized  equations  is  adequate  for  modeling  the 
physics;  or  subgrid  physics  cannot  be  ignored  and  local  mesh  refinement  must  be  used  to 
resolve  the  internal  structure  of  the  discontinuity.  There  are  problems  of  both  types,  and 
there  are  problems  with  multiple  length  scales  that  combine  these  types  (so  that  the 
discontinuity,  once  resolved,  contains  a  further  sharp  discontinuity  within  it). 
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Early  proposals  for  tracking  arc  described  in  Richtmyer  and  Morton  [20].  These 
ideas  were  further  developed  by  Moretti  [18].  These  solutions  are  of  high  quality  but 
appear  to  be  limited  to  somewhat  simplified  situations.  The  goal  of  the  work  reported 
here  is  to  implement  a  general  purpose  computational  package  based  on  front  tracking 
ideas  that  provides  highly  resolved  solutions  for  partial  differential  equations.  At  present 
it  is  applicable  to  problems  with  important  discontinuities  that  contain  ignorable  subgrid 
physics.  The  current  effort  is  based  on  systematic  use  of  data  structure  and  modular  pro- 
gramming concepts.  Our  implementation  is  designed  to  allow  flexible  adaptation  to  a 
variety  of  fluid  flow  problems.  Parallel  efforts  using  many  of  the  same  computational 
modules  have  been  directed  at  incompressible  two  phase  flow  in  oil  reservoirs  [10,7,8] 
and  instabilities  in  fluid  interfaces  [13]. 

In  this  paper  we  describe  the  algorithm  used  to  solve  a  hyperbolic  system  of  non- 
linear conservation  laws 

w,  +  Vf(w)  =  0.  (1.1) 

While  most  of  the  discussion  will  be  valid  for  general  systems  of  conservation  laws,  we 
will  be  concerned  specifically  with  the  Euler  equations  for  a  compressible,  invisdd, 
polytropic  gas,  for  which 


w  = 


m 


and  f(w)  = 


m 

mm 


+  P 


(1.2) 


Here  p  is  the  mass  density  of  the  fluid,  m  is  the  momentum  density  (m  =  pv  where  v  is 
the  fluid  velocity),  £  is  the  total  energy  density,  and  /?  is  the  thermodynamic  pressure, 
specified  by  the  polytropic  equation  of  state 

;,  =  (7-l)[£-|mp/2p]  (1.3) 

for  some  constant  "y  >  1.  These  equations  embody,  respectively,  the  conservation  of 
mass,  Newton's  law,  and  the  conservation  of  energy. 
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Wc  have  tested  the  rcsolutioii  of  the  front  tracking  niethod,  as  applied  to  gas  dynam- 
ics, by  solving  model  problems  on  coarse  or  moderate  grids.  The  results  of  these  tests  are 
described  in  Sec  6.  In  principle,  tracking  also  applies  to  complicated  problems  on  fine 
grids  if  fast  computers  are  used,  and  it  applies  to  problems  (such  as  combustion)  for  which 
subgrid  physics  is  not  ignorable  if  a  mesh  rcfmement  algorithm  is  incorporated,  but  these 
steps  have  not  yet  been  taken. 

2.   Concepts  and  Data  Stmctares 

In  this  section  we  describe  the  data  structures  used  m  our  implementation  of  front 
tracking.   The  discussion  follows  the  division  of  the  code  into  libraries. 

2.1.  Geometry  and  Interfaces 

The  basic  data  structure  defining  the  geoirsetry  and  topology  of  the  computational 
domain  is  called  an  interface.  An  interface  in  a  two-dimensional  domain  is  defined  as  a 
collection  of  oriented  curves.  A  carve  starts  and  ends  at  nodes,  and  consists  of  a  doubly 
linked  list  of  bonds  (so  that  each  bond  contains  a  pointer  to  the  next  and  the  previous 
bond).  Associated  with  each  node  and  with  the  beginning  and  end  of  each  bond  is  a 
point  that  defines  a  position  in  the  plane.  In  most  situations  the  curves  are  not  dlowed  to 
intersect  (except  at  nodes);  in  this  case  an  interface  divides  the  plane  into  connected  com- 
ponents. 

A  library  of  subroutines  is  devoted  to  supporting  these  data  structures  and  perform- 
ing elementary  operations,  including  creating  (i.e.  allocating,  storing,  and  installing), 
deleting,  copying,  printing,  and  reading  these  ob^xts.  Because  linked  lists  are  used,  the 
time  needed  to  add  or  delete  a  point  is  independent  of  the  number  of  points  on  the  inter- 
face. Other  routines  are  available  to  join  or  spHt  curves,  check  for  intersections,  zoom  on 
subdomains,  and  compute  the  topological  component  containing  a  given  point.  This  last 
computation  arises  frcquentiy  and  must  be  efficient.  By  suitable  preprocessing  of  the 
interface,  this  point  location  problem  can  be  made  effectively  independent  of  the  length  of 
the  interface  for  non-pathological  interfaces  [12]. 


The  interface  library  is  described  in  more  detail  in  a  separate  publication  [12]. 

2.2.  States  and  Fronts 

The  interface  data  structures  (and  the  programming  language  in  which  they  arc  writ- 
ten) allow  insertion  of  problem-specific  data  entries.  An  interface  equipped  with  data 
entries  specific  to  front  tracking  is  known  as  a  front.  A  front  contains  some  general  ideas 
of  physics  but  excludes  details  of  the  physical  equations  being  noodeled.  The  principal 
entries  specify  the  physical  nature  of  curves  and  nodes  and  the  state  variables  associated 
with  them. 

Curves  are  of  two  general  types:  boundary  curves  and  physical  curves.  A  boundary 
curve  is  further  classified  as  being  a  periodic  boundary,  a  Dirichlct  boundary  (coupling  to 
ambient  reservoir),  or  a  Neumann  boundary  (acting  as  a  reflecting  boundary).  The  vari- 
ous physical  types  of  curves  arc  defined  by  the  physics  being  modeled.  For  gas  dynamics 
the  physical  curves  arc  sound  waves  (including  shocks)  and  contact  discontinuities  (includ- 
ing material  discontinuities  and  slip  lines).  Nodes  arc  of  three  general  types:  fixed  nodes, 
boundary  nodes,  and  physical  nodes.  A  fixed  node  joins  boundary  curves;  for  example  it 
can  demarcate  a  comer  in  a  wall  or  separate  curves  corresponding  to  different  boundary 
conditions.  At  a  boundary  node  a  single  physical  curve  meets  a  boundary  curve,  such  as 
when  a  shock  moves  normally  along  a  wall.  The  classification  of  physical  nodes  for  gas 
dynamics  will  be  given  in  Sec  5. 

Throughout  most  of  the  code,  a  state  is  the  location  of  a  segment  of  storage  with  a 
fixed  size.  Exactly  what  information  is  stored  in  a  state  depends  on  the  physics  being 
modeled,  but  only  a  few  of  the  lowest  level  subroutines  (e.g.  the  Riemann  problem  solver 
and  the  Lax-Wendroff  elementary  integration  step)  depend  on  the  specific  form  of  this 
information.  For  gas  dynamics  there  are  four  real  numbers,  specifying  the  density,  two 
components  of  momentum  density,  and  the  energy  density,  together  with  the  address  of 
the  data  base  describing  the  equation  of  state.  Since  a  curve  in  the  front  represents  a  phy- 
sical wave  across  which  there  is,  in  general,  a  discontinuity  in  the  solution,  there  arc  two 
states  associated  with  each  point  of  a  curve,  corresponding  to  its  two  sides.  At  a  node 
there  is  a  state  associated  with  each  component  adjacent  to  it.  For  convenience,  however. 


we  associate  these  states  with  the  ends  of  the  curves  that  meet  at  the  node,  rather  than 
with  the  node  itself. 

The  front  library  is  a  collection  of  subroutines  that  support  states,  boundary  condi- 
tions, and  the  data  management  and  geometrical  aspects  of  the  dynamical  propagation  of 
curves  and  nodes.  All  physics-dependent  operations  are  confined  to  a  few  subroutines 
that  are  accessed  through  pointers  to  functions  that  arc  defined  in  a  physics  library. 
Briefly,  the  dynamical  propagation  of  the  front  is  as  follows.  For  each  curve  there  are 
two  sweeps  over  its  points:  the  first  employs  a  physics-dependent  subroutine  that  naovcs 
each  point  and  updates  its  associated  states  in  acconknce  with  waves  propagating  nor- 
mally to  the  curve;  the  second  constructs  a  three-point  stencil  of  states  centered  at  each 
point  and  calls  a  physics-dependent  subroutine  that  updates  the  center  state  by  taking 
account  of  waves  nKwing  tangentially  to  the  curve.  The  points  in  the  vidnity  of  each 
node  are  propagated  so  as  to  allow  for  interaction  of  the  waves  meeting  at  the  node.  For 
fixed  and  boundary  nodes  this  can  be  accomplished  using  the  subroutine  that  propagates  a 
point  together  with  geometrical  constructions;  for  physical  nodes  specialized  routines  arc 
called.  The  front  is  then  remeshed  according  to  one  of  several  algorithms  designed  to 
redistribute  points  to  enhance  curvature  resolution  and  to  prevent  pile-ups  and  thinning. 
To  define  the  states  at  new  front  mesh  points  the  rcmeshing  subroutine  uses  a  physics- 
dependent  subroutine  to  interpolate  between  nearby  states  belonging  to  the  same  com- 
ponent. Fmally  the  front  is  checked  for  self-intersections.  If  intersections  are  found,  a 
physics-dependent  subroutine  that  untangles  the  front  is  called.  Further  details  on  the 
propagation  of  the  front  will  be  given  in  Sec  4. 

2.3.  Interior  Waves 

The  region  away  from  the  curves  defining  the  front  is  called  the  interior.  Since  the 
curves  represent  possible  discontinuities  in  the  solution,  each  component  of  the  interior  is 
regarded  as  the  domain  for  a  separate  initial/boundary-value  problem.  Depending  on  the 
physics  being  modeled,  elliptic,  parabolic,  or  hyperbolic  equations  are  to  be  solved.  For 
gas  dynamics  the  system  is  hyperbolic,  so  we  will  restrict  our  attention  here  to  a  library  of 
routines  for  solving  hyperbolic  equations,  especially  hyperbolic  conservation  laws. 


State  variables  in  the  interior  are  associated  with  points  on  a  rectangular  grid.  Each 
grid  point  belongs  to  a  specific  component.  Tbe  value  of  the  solution  at  an  arbitrary  posi- 
tion in  the  computational  domain  is  obtained  through  bilinear  interpolation  of  interior 
states  or,  in  case  the  position  is  close  to  the  front,  through  linear  interpolation  between 
interior  states  and  states  on  the  front  belonging  to  the  same  component.  The  calculation 
of  such  states  is  necessary  in  the  normzd  propagation  sweep  along  the  front,  so  care  is 
taken  to  divide  each  component  into  triangular  elements  that  enable  efficient  interpola- 
tion. 

To  update  the  interior  states  in  time  a  number  of  schemes  arc  provided.  At  present 
the  random  choice,  scalar  upwind,  and  Lax-Wcndroff  schemes  using  spatial  operator  split- 
ting, as  well  as  a  fully  two-dimensional  Lax-Wendroff  scheme,  have  been  implemented. 
It  is  critical  to  the  front  tracking  method  that  the  states  used  in  the  stencils  for  these 
schemes  all  belong  to  the  same  component.  In  our  application  to  gas  dynamics  we  have 
used  the  Lax-Wendroff  schemes,  which  we  describe  further  in  Sec.  3. 

2.4.  Utilities  and  Drivers 

Various  general  purpose  routines  for  storage  allocation,  debugging,  and  input/output 
are  contained  in  a  utility  library.  Routines  that  govern  the  control  flow  of  the  time  loop, 
initialization,  diagnostic  analysis  of  the  solution,  and  printing  are  contained  in  a  driver 
library.  Also  in  the  driver  library  are  routines  that  specially  format  the  input  and  output 
to  interface  to  a  package  of  data  analysis  and  graphics  programs. 

2.5.  Physics  Snbrontines 

As  indicated  above,  certain  operations  on  state  variables  and  points  on  the  front  ulti- 
mately depend  on  the  physics  being  modeled.  All  of  these  features  have  been  isolated  in 
a  handful  of  subroutines  that  are  accessed  though  abstract  function  pointers.  For  gas 
dynamics  the  key  subroutines  needed  are:  the  subroutine  that  propagates  a  point,  which  is 
essentially  a  non-local  Riemann  problem  solver;  the  subroutine  that  updates  the  states  dur- 
ing the  tangential  sweep  of  the  front,  which  is  a  one-dimensional  L.ax-Wendroff  operator; 
the  specialized  subroutines  that  propagate  physical  nodes,  which  solve  two-dimensional 
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Riemann  problems;  and  the  stencil  operator  for  the  interior  scheme,  which  is  a  one-  or 
two-dimensional  Lax-Wendroff  operator.  These  subroutines  will  be  described  in  the  fol- 
lowing. In  addition  there  is  the  subroutine  that  initializes  the  front  and  interior  according 
to  the  test  problem  being  studied  along  with  miscellaneous  subroutines  for  diagnostic 
analysis  and  printing  of  the  solution. 

3.  The  Interior  Scheme 

The  numerical  approximation  to  the  solution  of  the  non-linear  conservation  law  (1.1) 
is  represented  by  state  variables  at  the  nodes  of  a  rectangtilar  grid  together  with  double- 
valued  state  variables  at  the  points  on  an  interface.  The  solution  at  an  arbitrary  position 
in  the  computational  domsun  is  obtained  by  interpolating  among  these  interior  and  front 
states,  as  we  now  describe.  The  interpolation  position  lies  in  some  rectangular  grid  block 
defined  by  four  grid  nodes;  it  is  also  regarded  as  lying  in  some  specified  topological  com- 
ponent. Since  the  front  represents  a  possible  discontinuity  in  the  solution,  it  is  important 
not  to  interpolate  between  states  corresponding  to  different  components.  K  the  four  grid 
nodes  all  lie  in  the  same  component  as  specified,  then  we  may  obtain  the  solution  using 
bilinear  interpolation.  If,  however,  not  all  of  the  grid  nodes  lie  in  the  same  component, 
the  front  cuts  through  the  grid  block  to  separate  the  grid  nodes  from  one  another.  In  this 
case  only  the  grid  nodes  that  lie  in  the  specified  component  of  the  interpolation  position, 
togetlKr  with  states  on  the  appropriate  side  of  the  front,  are  used  in  the  interpolation. 
Thus  we  triangulate  each  grid  block  that  is  crossed  by  the  interface  and  use  linear  interpo- 
lation on  these  triangvilar  elements.  Our  present  triangulation  scheme  determines  the 
points  on  the  front  where  it  crosses  the  vertical  and  horizontal  grid  lines;  these  points, 
along  with  the  grid  nodes  lying  in  the  correct  component,  define  the  comers  of  the  trian- 
gles. In  this  scheme  it  is  computationally  efficient  to  determine  the  triangle  in  which  the 
interpolation  point  lies,  since  in  all  but  one  case  the  triangles  are  zirranged  in  a  star-like 
fashion  about  one  common  comer. 

For  the  calculation  of  the  solution  in  the  interior  regions  we  use  either  an  operator 
split  one-dimensional  Lax-Wendroff  or  a  fully  two-dimensional  Lax-Wendroff  scheme. 
Here  we  discuss  only  the  two-dimensional  (imsplit)  scheme.   The  two-dimensional  scheme 
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involves  two  half  steps.  To  facilitate  the  coupJu-J  of  i'l^  front  and  interior,  the  front  is 
also  advanced  in  half  steps.  Since  the  Lax-Wendroif  schciiic  is  a  leapfrog  composition  of 
two  Lax-Friedrichs  steps  it  is  enough  to  describe  the  Lax-Friedrichs  scheme.  This  scheme 
usually  assumes  the  initial  data  to  be  known  at  the  four  comers  of  a  square.  In  our  appli- 
cation, however,  there  are  irregular  squares,  i.e.  ones  for  which  one  or  more  of  its  comers 
is  cut  off  by  the  front  and  thus  lies  in  the  wrong  component.  To  circumvent  this  difficulty 
we  view  the  Las-Friedrichs  scheme  as  defined  by  a  flux  balance:  the  sum  of  the  fluxes 
through  the  sides  of  a  mesh  block  determines  the  change  in  time  of  a  conserved  quantity 
integrated  over  the  block.  From  this  point  of  view  the  Lax-Friedrichs  scheme  is  dcfmed 
for  irregular  squares  as  well  as  regular  squares.  In  fact,  the  propagation  of  the  front  also 
determines  the  fluxes  through  the  front  and  thereby  through  the  sides  of  the  irregular 
squares.  The  triangulation  that  is  constructed  for  interpolating  the  solution  is  used  to 
define  these  irregular  squares. 

To  be  more  explicit,  consider  a  particular  grid  block  with  side  length  Ax.  Given  the 
solution  at  time  /q.  we  wish  to  calculate  the  solution  at  the  center  x^  o^  this  block  at  time 
ti  =  Iq  +  Ar/2.  The  center  of  the  grid  block  lies  in  some  topological  component  at  the 
later  time  t^,  and  it  is  the  solution  in  this  component  that  we  wish  to  evaluate.  Let  A(/) 
denote  the  portion  at  time  /  of  the  grid  square  lying  in  this  component,  and  let  \A(t)\ 
denote  its  area.  Then  by  integrating  the  conservation  law  (1.1)  over  this  region  and  then 
integrating  in  time  from  IqIo  t-^  we  obtain 

f  w(x,/i)  =    /  w(x,/o)  +  Jdt  f  a(x,/)w(x,0  -  fdt  f  /i(x,/)f(w(x,0).  (2.1) 

M',)  M'o)  'a      3A{t)  f„      dA{t) 

Here  n(x,t)  and  a(x,t),  for  points  x  on  the  boundary  dA(t),  are  the  unit  normal  vector 
and  the  the  normal  speed  of  the  boundary.  Thus  for  points  x  on  the  tracked  ircnt,  ct(x,/) 
is  the  normal  speed  of  the  front,  while  for  other  points  it  vanishes.  In  order  to  obtain  a 
numerical  approximation  to  this  formula  that  is  stable,  we  note  the  identity 

Jdt  f  a(3,0w  =  (  lA(/i)|  -  !.l(/o)!  )  w  (2.2) 

for  any  constant  state  w.   By  suitably  choosing  w  we  can  arrange  that  the  integral  on  the 
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left  in  (2.2)  approximates  the  secxmd  integral  un  ibc  right  in  (2.1)  with  relative  error 
O(Ax^).  For  instance,  if  the  front  at  time  t^  consists  of  a  single  line  segment  crossing  the 
grid  block,  then  we  can  take  w  to  be  the  solution  on  the  front  at  the  center  of  this  seg- 
ment.  By  approximating  the  other  integrals  in  (2.1)  we  obtain  the  formula 

lA(/i)|w(x„/i)  =  kl(/o)|wo  +  (  \A{h)\  -  \A{to)\  )  w 

-  "T    /  n(x,ro)f(w(x,/o))  +  OiM^Ax")  +  0(A/Ax-^)    (2.3) 

^   aA(g 

in  which  wq  represents  the  average  value  of  the  solution  over  the  region  A(/o),  so  that  it 
corresponds  to  the  correct  component.  If  \A(ti)\  /Ajtp-  is  not  too  small,  equation  (2.3)  pro- 
vides a  suitable  approximation  for  w(Xf  ,/i);  otherwise  x^,  is  in  general  so  close  to  the  front 
that  it  is  appropriate  to  approximate  w(x<,,/i)  with  states  on  the  front.  We  note  that  when 
the  front  does  not  cross  the  grid  square,  equation  (2.3)  reduces  to  the  usual  Lax- 
Friedrichs  approximation.  Finally,  the  solution  at  time  tQ  +  ^t  is  obtained  using  two 
Lax-Friedrichs  steps  in  a  leapfrog  combination. 

4.  The  Front  Scheme 

The  propagation  of  the  front  involves  the  motion  of  points  on  the  front  and  the  evo- 
lution of  the  states  on  the  front.  The  propagation  of  points  in  the  vicinity  of  nodes,  where 
in  general  many  curves  meet,  is  described  in  Sec.  5;  here  we  consider  points  in  the  interior 
of  a  curve. 

Let  ZQhc  an  interior  point  on  the  front  at  time  Iq.  Points  z^^^  =  ^^'^\^)  near  zq  are 
labeled  by  their  arclength  displacement  s  away  from  zg  =  z(°)(0).  A  coordinate  system 
for  a  neighborhood  of  Zq  is  constructed  as  follows  (sec  Fig.  1).  Each  pair  (r,  s) 
corresponds  to  the  point  x(r,  s)  that  is  a  distance  r  along  the  line  drawn  normal  to  the 
front  through  the  point  z^°^(^).  The  coordinate  curves  s  =  const,  are  thus  lines  perpendic- 
ular to  the  front.  Let  n(r,  s)  and  s(r,  j)  be  the  parallel  translates  to  x(r,  s)  of  the  unit 
vectors  normal  and  tangent  to  the  front  at  z^°)(j),  respectively.   The  system 

w,  +  Vf(w)  =  0  (4.1) 
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of  conservation  laws  may  be  written  in  terms  of  the  normal  and  tangential  derivatives  /i-V 
and  i-V  as 

w,  +  n[  (nV)  f(w)  ]  +  s[  (iV)  f(w)  ]  =  0.  (4.2) 

We  wish  to  solve  this  system  subject  to  the  the  initial  conditions  w(/  =  Iq)  =  \r^^\ 
where  w*^*^^  is  smooth  except  for  possible  jump  discontinuities  across  z^^^  Let  w^'^  denote 
the  exact  solution  at  time  Iq  +  At  to  this  problem;  w^''  is  smooth  except  for  possible 
jumps  at  the  exact  front  position  z^''  =  z^'\s).  We  say  that  a  computed  solution  w^"^^  is 
correct  at  the  approximate  front  z^*^)  =  z^'^\s)  through  order  At  if 

z<^%)  =  z^%)  +  0{At^)  (4.3) 

and 

w(')(z(')(-y),  tQ  +  At)  =  ^r^'\z<^%),  tQ  +  At)  +  0(At^),  (4.4) 

where  by  this  last  equation  we  mean  equality  of  the  states  on  corresponding  sides  of  the 
jump  discontinuities. 

The  computed  solution  of  system  (4.1)  is  obtained  using  operator  splitting  of  the  sys- 
tem (4.2)  in  the  normal  and  tangential  directions.   First  the  normal  equations 

w,  +  n-[  (/i-V)  f(w)  ]  =  0,   w(/  =  to)  =  w^O)  (4.5) 

are  solved  approximately  to  obtain  the  computed  front  position  z^*^^  =  z^'^\s)  and  the  com- 
puted normal  solution  w^'^)  at  time  /q  +  ^f-  Then  the  tangential  equations 

w,  +  s[  (iV)  f(w)  ]  =  0,  w(/  =  tQ)  =  w(^)  (4.6) 

are  solved  approximately  to  obtain  the  solution  w^'^^  at  time  /q  +  A/. 

Let  us  analyze  the  error  involved  in  this  procedure.  The  normal  front  speed  a(s)  at 
time  tQ  satisfies 

z^%)  =  z^%)  +  At  <t(s)  n(0,  s)  +  0(A/2).  (4.7) 

By  Eq.  (4.2)  and  Taylor's  theorem, 

w(')(z(')(5).  tQ  +  At)  =  w(0)(z(0)(^),  to)  +  At  fjis)  [  (n-V)  w(0)  ]{z^%),  to) 
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-  A/  n[  (nV)  f(w(0')  ](z(%),  /q) 

-  ^t  s[  (i-V)  f(w(0))  ](z(0)(^),  /o)  +  0(A/2),  (4.8) 

Assume  that  the  computed  normal  solution  w^*^^  is  a  correct  solution  at  the  approximate 
front  z^'^^  =  T,^'^\s)  to  the  normal  equations  (4.5)  through  order  A/.  Since  the  front  posi- 
tion obtained  by  solving  the  normal  equations  will  agree  through  order  At  with  the  exact 
front  position,  we  conclude  that 

w^^'(z(^'(^),  /o  +  AO  =  w(0)(z(0)(^).  to)  +  A/  u{s)  [  (n-V)  w^O)  ](z^%),  /q) 

-  A/  rt[  (/iV)  f(w^°))  ](z(°)(5),  /o)  +  0(A/2).  (4.9) 

Assume  in  addition  that  w^"^^  is  a  correct  solution  at  the  approximate  front  z^'^^  =  z^'^\s) 
to  the  tangential  equations  (4.6)  through  order  A/.   Therefore 

w^'^)(z(^)(^),/0  +  AO  =  w(^)(z(^)(5),  tQ  +  At) 

-  At  s[  (i(j)-V)  f(w(^))  ]iz^%),to  +  At)  +  0(At^).  (4.10) 

Combining  these  last  two  equations  shows  then  that  w^*^'  is  a  correct  solution  at  the 
approximate  front  z^*^^  =  z^'^\s)  to  the  full  equations  (4.2)  through  order  At.  In  the  fol- 
lowing we  demonstrate  how  normal  and  tangential  solutions  satisfying  the  above  assump- 
tions can  be  obtained. 

First  we  propagate  the  front  and  the  states  along  the  front  in  the  normal  direction. 
The  solution  at  time  /q  is  evaluated  on  both  sides  of  the  front  at  zq,  yielding  a  left  and  a 
right  state.  The  solution  at  time  /q  is  also  evaluated  at  a  normal  distance  Ar  on  each  side 
of  the  front.  These  states  will  be  used  to  calculate  the  waves  that  impinge  on  the  front 
from  the  interior.  Notice  that  if  the  front  contains  curves  too  dose  to  each  other,  these 
new  evaluation  points  zg  ±  Ar-n  may  be  in  different  components  from  the  points  zg  ±  0-n 
at  the  front.  In  this  case  the  evaluation  point  is  shifted  into  the  correct  component. 
(Conceptually,  the  state  at  a  point  outside  a  given  component  is  obtained  by  extrapola- 
tion.) Thus  we  always  have  as  data  two  left  states  corresponding  to  one  component  and 
two  right  states  corresponding  to  a  second  component. 
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These  states  are  used  as  initial  data  for  an  extended  or  non-local  Riemann  problem. 
Higher  order  solutions  of  these  Riemann  problems  have  been  discussed  before.  To  reduce 
sampling  error  in  the  random  choice  method,  a  steady  state  Ansatz  has  been  used  [6,11] 
to  extend  local  Riemann  data  over  a  mesh  block,  thereby  obtaining  a  non-local  Riemann 
problem,  which  is  solved  to  higher  order.  Higher  order  Godunov  schemes  [1,2]  also 
employ  ideas  related  to  solutions  of  non-local  Riemann  problems.  In  the  present  scheme 
we  solve  the  non-local  Riemann  problem  as  follows. 

Using  the  left  and  right  states  located  at  the  front  we  solve  an  ordinary  Riemann 
problem.  The  solution  is  the  correct  answer  to  the  non-local  problem  at  time  /q  +  0.  and 
it  is  used  to  approximate  the  propagation  speeds  of  the  characteristics  backward  from 
/q  -I-  A/  to  tQ.   We  find  starting  points  for  these  characteristics  in  the  normal  intervals 

[  zo  ±  0-n,  Zq  ±  Ar-n  ].  (4.11) 

The  states  at  these  points  are  calculated  using  linear  interpolation  between  the  states  en 
the  front  and  the  states  at  the  normally  displaced  points.  In  this  way  we  determine  which 
waves  from  the  normal  intervals  enter  the  front.  Using  differential  equations  in  charac- 
teristic form  together  with  the  Rankine-Hugoniot  jump  conditions  we  compute  the  states 
to  be  associated  with  the  propagated  front  at  time  /g  +  A/  . 

The  solution  of  the  normal  sweep  along  the  front  is  taken  as  initial  values  for  the 
tangential  sweep.  By  linear  interpolation,  the  states  on  each  side  of  the  front  can  be 
defined  everywhere  along  the  front.  Therefore  we  can  evaluate  the  normally  propagated 
solution  at  each  mesh  point  and  at  two  neighboring  points  displaced  a  distance  Aj  along 
the  front.  Using  these  three  stencil  points  and  the  one-dimensional  Lax-Wendroff 
scheme,  adapted  for  the  tangential  equations  (4.6),  we  determine  the  tangentiaUy  pro- 
pagated state  variables.  Notice  that  tangential  propagation  of  the  points  on  the  front  is 
equivalent  to  a  remeshing  of  the  front,  in  the  limit  A/  -  0,  so  it  is  not  essential  to  move 
these  points  during  the  tangential  sweep. 

In  what  follows  the  normal  sweep  is  described  in  more  detail  for  gas  dynamics. 
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4.1.  The  Normal  Sweep  for  Shock  Waves 

There  are  two  types  of  shocks  for  the  gas  dynamic  equations:  backward  and  forvi'ard 
(or  1  and  3)  shocks.   Here  a  t -shock  is  a  discontinuity  satisfying 

Xi_i(w,)   <    a   <    Xi(vr;)     and    X^Cw,)   <    a   <    Xi  +  i(w^),  (4.12) 

where  a  "<;  the  normal  propagation  speed  of  the  shock,  and  X^  =  u  —  c,  X2  =  m,  and 
X3  =  u  +  r  are  the  eigenspeeds  of  normal  equations  (4.5),  u  =  n-m/p  being  the  normal 

velocity  and  c  =  {"ipfpY  being  the  speed  of  sound  [21].  Thus  there  are  three  characteris- 
tics entering  each  point  on  the  left  (rig'  ^  side  and  one  characteristic  entering  each  point 
on  the  right  (left)  side  of  a  backward  (foiward)  shock.  We  therefore  solve  three  charac- 
teristic equations  to  obtain  the  state  ahead  of  the  shock,  and  then  solve  one  characteristic 
equation  together  with  the  Ran^ine-Hugoniot  conditions  to  obtain  the  state  behind  the 
shock.   We  shall  describe  the  algorithm  for  a  forward  shock  only. 

The  normal  equations  (4.5)  for  gas  dynamics,  as  written  in  characteristic  form,  are 

d   f   2c     _       _   c    d 


and 


d   ,_2£_        N  _   <^    d 


d\-}    y  —  l  7  d\2 

where  V  =  s-m/p  is  the  tangential  velocity  and  S  =  log(p/p'^/{y-l)  is  the  entropy. 
There  are  no  waves  transmitted  to  the  right  side  of  a  forward  shock,  so  the  characteristic 
equations  determine  the  state  on  the  right  side.  We  obtain  an  approxinjate  solution  of  the 
characteristic  equations  by  solving  the  following  difference  equations: 

2cr  2ci  {Cf  +  c{) 

+  "1  =  ^7 {^r  -  -S^l). 


-y-1  '■        -y-1  ^  2-y 

S,  =  52,    V,  =  V2,  (4.14) 
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where  the  subscript  r  refers  to  values  at  time  /g  +  ^^  on  the  right  side  of  the  shock,  and 
the  subscripts  1,  2  and  3  refer  to  values  at  the  feet  of  the  X^,  X2,  and  X3  characteristics  at 
time  /q  (see  Fig.  2). 

To  obtain  these  states  at  the  feet  of  the  characteristics  we  first  solve  the  Riemann 
problem  between  the  left  and  right  states  on  the  front  at  time  /q  to  obtain  a  preliminary 
shock  position  at  time  /q  +  ^t-  Then  the  three  characteristics,  approximated  as  straight 
lines,  are  drawn  backward  from  the  shock  position  at  /g  +  ^^  to  foot  positions  in  the  nor- 
mal interval  [zq,iq  +  Ar-n  ]  at  time  Iq.  The  corresponding  states  are  obtained  using 
linear  interpolation  between  the  right  state  at  zq  and  the  state  at  zq  +  Ar-ii.  Using  these 
states  in  the  above  difference  approximation  to  the  characteristic  equations  yields  a  right 
state  at  time  Iq  +  A/  that  is  correct  through  order  At. 

On  the  left  side  of  the  shock  only  the  X3  characteristic  impinges  on  the  shock.  The 
corresponding  characteristic  equation  is  approximated  by  the  difference  equation 
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where  the  subscript  3/  refers  to  values  at  the  foot  of  the  left  X3  characteristic.    The 
Rankine-Hugoniot  conditions  for  a  forward  shock  are 
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Given  the  solution  obtained  above  for  the  state  W;.  on  the  right  side  of  the  shock,  the 
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characteristic  and  Rankine-Hugoniot  equations  yield  an  approxknate  left  state  at  time 
Iq  +  A/  that  is  correct  through  order  A/.  Finally,  the  propagation  speed  of  the  shock  is 
calculated  by  averaging  the  shock  speeds  at  time  Iq  and  at  time  /q  +  ^f- 

4.2.  The  Normal  Sweep  for  Contact  Discontinuities 

A  contact  discontinuity  separates  two  states  having  different  values  of  density  and 
tangential  velocity  but  the  same  values  of  pressure  and  normal  velocity.  The  propagation 
speed  of  a  contact  is  the  normal  particle  speed  u ,  so  only  one  family  of  characteristics  on 
each  side  enters  the  discontinuity.   Thus  we  have 

d    ,   2c  -        c     d    „         .       d    .   2c      ,     .        c     d    _  /.^-n 

— -( -  u)  =  --^z-S    and    -—( +  u)  =  —7—5  (4.17) 
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on  the  right  and  left  sides  of  the  contact,  respectively.  In  addition  the  entropy  and  the 
tangential  velocity  on  each  side  of  the  contact  is  constant. 

The  corresponding  difference  equations  are 
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where  the  subscripts  /  and  r  refer  to  left  and  right  values,  respectively,  at  time  /q  +  ^t, 
while  the  subscripts  Ir  and  2r  refer  to  values  at  the  feet  of  the  right  \i  and  X2  characteris- 
tics at  time  /q,  and  the  subscripts  2/  and  3/  refer  to  values  at  the  feet  of  the  left  X2  and  X3 
characteristics  at  time  /q  (see  Fig.  3).  These  states  at  the  feet  of  the  characteristics  arc 
obtained  by  linear  interpolation  in  the  same  fashion  as  described  for  shock  waves.  In 
addition  to  the  characteristic  equations  we  have  the  Rankine-Hugoniot  conditions  that 
u/  =  u^  and  pi  =  Pr  for  a  contact.  By  solving  these  equations  we  obtain  an  approximate 
solution  to  the  normal  equations  that  is  correct  at  the  front  through  order  A/.  The  propa- 
gation speed  of  the  contact  is  calculated,  as  for  shocks,  by  averaging  the  speeds  at  time  /q 
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and  at  time  Iq  +  A/. 

4.3.   The  Normal  Sweep  for  Boundaries 

The  boundaries  are  considered  as  part  of  the  front,  and  the  states  on  the  boundaries 
are  updated  using  the  normal/tangential  operator  splitting  method.  The  boundary  condi- 
tions under  consideration  are  either  Dirichlet  or  Neumann.  Periodic  boundary  conditions 
are  implemented  so  as  to  have  no  effect  on  the  solution  and  are  not  discussed  further. 

Dirichlet  boundary  conditions  for  the  gas  dynamics  equations  are  defined  physically 
as  the  coupling  to  an  ambient  reservoir  (inlet  or  outlet).  Mathematically,  they  are  cfefined 
by  the  specification  of  a  boundary  state.  This  state  plays  the  role  of  far-field  conditions. 
Information  about  the  far-field  conditions  are  propagated  into  the  computational  region 
only  by  incoming  characteristics. 

Let  w/  and  w^  be  the  left  and  right  states  at  a  given  point  zq  on  the  boundary.  Let 
w//  be  the  state  at  a  normal  distance  Ar  on  the  left  side,  and  w^  the  state  at  a  normal  dis- 
tance Ar  on  the  right  side  of  the  boundary.  If  the  exterior  of  the  computational  region  is 
on  the  left  (right),  w/;  (w^)  is  the  far-field  boundary  condition.  To  obtain  the  normally 
propagated  solution  at  the  boimdary  we  first  perform  a  Lax-Wendroff  step  using  w;/,  w/, 
and  w/  as  the  values  at  the  points  Zq  —  ^r-n,  Zq,  and  Zq  +  Ar-n,  respectively,  obtaining 
the  intermediate  state  w|'^  We  do  the  same  on  the  right  side  of  the  curve  obtaining  w^'^ 
Then  the  computed  normal  left  and  right  states  w/'^'  and  w^'^^  are  obtained  by  solving  the 
Riemann  problem  between  the  intermediate  states  wf'^  and  w;'^. 

At  Neumann  boundary  curves  we  require  that  the  normal  velocity  vanish.  Thus  only 
one  sonic  characteristic,  X^  or  X3,  impinges  on  the  wall.  Suppose  the  computational 
region  is  on  the  right  side  of  the  boundary.  Just  as  in  the  case  of  a  contact  discontinuity, 
but  taking  into  account  that  the  normal  velocity  u  vanishes  at  the  boundary,  the  charac- 
teristic equations  may  be  approximated  by  the  difference  equations 

5  =  ^2,    V  =  V2, 

and 
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+  "1  =  -^(-^  -  "^l).  (4-19) 


-y-1         7-1  '  2-y 

where  the  subscripts  1  and  2  refer  to  values  at  the  feet  of  the  X^  and  X2  characteristics  at 
time  /q.  These  equations  may  be  solved  to  yield  a  solution  to  the  nonnal  equations  that  is 
correct  through  order  A/. 

5.  Two-Dimensional  Riemann  Problems 

An  intersection  of  two  or  more  curves  of  an  interface  is  called  a  node.  In  the  con- 
text of  front  tracking  a  node  is  the  center  of  the  interaction  of  the  waves  meeting  at  a 
point.  In  contrast  to  the  problems  discussed  in  Sec.  4,  the  geometry  near  a  node  does  not 
necessarily  admit  resolution  into  normal  and  tangential  components,  so  its  evolution  can- 
not be  determined  by  solving  a  sequence  of  one-dimensional  problems. 

A  node  is  categorized  by  the  (circular)  ordering  and  the  types  of  the  waves  emanat- 
ing from  it.  In  the  near  vicinity  of  a  node  the  curves  are  approximately  straight  lines  that 
separate  wedge-shaped  regions.  Thus  we  are  led  to  define  a  two-dimensional  Riemann 
problem  to  be  an  initial  value  problem  for  a  two-dimensional  conservation  law  having  data 
that  is  either  piecewisc  constant  or  is  a  simple  centered  rarefaction  wave  in  each  of  a  finite 
number  of  wedges.   The  form  of  the  data  is  thus 

w  =  w^      for  Oy-i  <  e  <  Oy  (5.1) 

for  y=l,...,n,  where  Qq=^„,  and  each  wy  is  either  a  constant  or  a  centered  wave.  Such 
problems  have  been  studied  for  scalar  conservation  laws  [14,16,23],  but  only  special  solu- 
tions are  known  for  systems  of  conservation  laws.  In  analogy  with  the  solution  of  one- 
dimensional  Riemann  problems,  the  solution  of  a  two-dimensional  Riemann  problem  will 
evolve  into  a  more  complicated  configuration  containing  several  nodes  (two  dimensional 
elementary  waves  interactions)  moving  apart  from  one  another. 

We  simplify  the  problem  in  two  ways.  First  we  only  look  for  elementary  wave 
interactions.  These  are  recognized  as  being  stable  under  forward  time  evolution  and  are 
associated  with  a  single  node  in  the  solution  of  a  general  Riemann  problem.  Secondly  we 
restrict  ourselves  to  problems  that  are  generic  under  change  of  initial  conditions,  and  to 
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solutions  that  are  limits  of  viscous  solurions  as  the  viscosity  parameter  tends  to  zero.  At 
this  point  physical  intuition  and  a  knowledge  of  experimental  facts  allows  a  guess  concern- 
ing a  (nearly  complete?)  list  of  elementary  wave  interactions.  Even  under  further  restric- 
tions (such  as  a  pwlytropic  gas  with  only  small  amplitucfe  waves)  there  seems  to  be  no 
mathematical  analysis  of  the  problem  of  classifying  elementary  waves  in  two-dimensional 
Riemann  problems. 

As  mentioned  in  Sec.  2.2,  a  node  is  either  a  fixed  node,  a  boundary  no<fc;,  or  a  phy- 
sical node.  At  a  fixed  node  only  bouncbry  curves  ixiect.  Here  we  describe  the  case  when 
only  two  boundary  curves  meet.  Such  a  node  is,  for  exzimple,  a  comer  of  a  wall  (when 
both  curves  are  of  Neumann  type)  or  the  end  of  an  inlet  or  outlet  (where  there  is  a  transi- 
tion from  Neumann  to  Dtrichlet  boundary  conditions).  Its  propagation  involves  only  the 
dynamics  of  the  state  variables  associated  with  the  node.  In  the  case  of  orthogonal  curves 
the  updated  states  arc  determined  using  operator  splitting  in  the  directions  normal  and 
tangential  to  one  of  the  curves.  For  transitions  from  Neumann  to  Dirichlet  boundaries 
meeting  at  general  angles  this  algorithm  should  also  be  correct.  In  the  case  of  an  arbi- 
trary angle  between  walls  we  use  a  scheme  in  which  flux  through  a  polygonal  element 
near  the  comer  is  balanced  to  obtain  the  states  around  the  node. 

The  general  phenomena  at  a  node  where  a  physical  curve  interacts  with  a  boundary 
curve  appears  to  be  as  follows.  If  the  boundary  curve  is  of  Dirichlet  or  periodic  type,  no 
interaction  between  Lhe  boundary  and  the  physical  wave  occurs.  Thus  normal/tangential 
splitting  is  correct  so  long  as  the  states  at  positions  outside  of  the  computational  region  are 
obtained  by  applying  the  appropriate  boimdary  condition.  We  therefore  restrict  the  dis- 
cussion to  the  interaction  of  a  physical  wave  with  a  Neumann  boundary,  i.e.  a  wall.  For 
a  general  conservation  law,  any  physical  wave  can  meet  a  fi-t  wall  normally;  again 
normal/tangential  spUtting  is  correct.  For  gas  dynamics  a  single  wave  meeting  a  wall 
obliquely  can  only  be  a  contact  discontinuity:  a  slip  line  can  meet  a  wall  only  tangcntially, 
and  a  shock  wave  can  only  meet  a  wall  normally,  since  the  normal  velocity  at  the  wall 
must  be  zero.  Configurations  in  which  two  or  more  waves  meet  at  the  wall  can  occur, 
however.  In  the  regular  reflection  of  a  shock  wave,  two  shock  waves  meet  at  a  wall,  one 
being  incident  (incoming  to  the  wall),  the  other  being  reflected.   A  shock  can  also  meet  a 
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smooth  portion  of  a  wall  obliquely  in  the  presence  of  a  slip  line.  In  effect  the  slip  is 
equivalent  to  a  comer  in  the  wall  and  the  oblique  shock  is  an  attached  bow  shock.  Such 
flow  can  occur  in  supersonic  nozzle  flow  and  leads  to  boundary  separation.  Additional 
phenomena,  for  instance  the  attached  bow  shock,  occur  in  the  presence  of  discontinuities 
in  the  boundary,  such  as  at  a  comer  of  a  wall  or  at  an  inlet.  In  the  case  of  one  physical 
wave  the  propagation  of  the  node  is  obtained  by  operator  splitting  in  the  normal  and 
tangential  dinections  to  the  physical  curve.  In  the  case  of  two  waves,  one  incident  and  the 
other  reflected,  we  use  the  Rankine-Hugoniot  conditions  for  the  incident  wave  and  shock 
polar  calculations  [3]  to  obtain  the  angle  of  reflection  and  the  states  around  the  node.  As 
with  the  one-dimensional  Riemann  problem,  non-local  information  should  be  (and,  in 
nK)st  cases  implemented,  is)  used  to  compute  the  effects  of  waves  entering  or  transmitted 
through  the  node. 

The  same  ideas  apply  to  the  classification  of  interior  nodes.  In  fact  a  contact  discon- 
tinuity serves  as  a  reflection  surface  for  shock  waves,  so  that  one  or  two  shock  waves  may 
appear  on  each  side  of  the  contact.  The  sound  speeds  on  either  side  of  the  contact  in  gen- 
eral differ.  Moreover  the  speed  of  an  incident  shock  as  it  moves  along  a  contact  depends 
on  the  angle  it  makes  with  the  contact.  Using  both  sound  speeds  and  incident  angles,  a 
fast  side  and  a  slow  side  of  the  contact  can  be  identified.  On  the  fast  side  the  wave  confi- 
guration is  either  a  normal  shock  or  a  regular  reflection,  while  on  the  slow  side  a  single 
transmitted  shock  occurs.  Because  the  contact  can  bend  at  the  point  of  incidence  of  the 
shock,  the  reflected  wave  can  be  a  rarefaction  [5].  Also,  two  incident  shock  waves  can 
meet  each  other  producing  in  addition  a  pair  of  reflected  waves  and  in  general  a  trailing 
contact  discontinuity.  However  if  the  two  waves  belong  to  the  same  family,  their  meeting 
provides  a  transmitted  shock,  a  contact,  imd  a  reflected  rarefaction  wave.  The  Mach  stem 
configuration  is  yet  another  node  type.  A  three  wave  interaction  is  the  generic  configura- 
tion for  contact  waves.  The  classification  becomes  more  complicated  when  waves  ending 
at  points  of  zero  strength  are  allowed.  The  propagation  of  these  interior  nodes  is  not 
implemented  at  present. 
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6.  Validatioii 


We  present  four  series  of  calculations  that  test  the  two-dimensional  propagation  of 
discontinuous  waves.  The  modes  of  wave  propagation  in  two-dimensional  gas  dynamics 
are  shock  waves  (nonlinear  sound  waves),  contact  discontinuities  (temperature  fronts), 
and  slip  lines  (transverse  velocity  discontinuities).  The  problems  testing  these  modes  are 
chosen  to  allow  verification  by  means  of  indepencfent  solutions. 

6.1.   Circularly  Symmetric  Problems 

In  our  first  series  of  tests  we  study  the  propagation  of  curved  shock  waves  Jind  con- 
tact discontinuities  evolving  from  circularly  symmetric  initial  data.  An  initial  pressure 
discontinuity  across  a  circle  gives  rise  to  two  discontinuous  waves,  a  shock  and  a  contact. 
For  these  problems  the  comparison  solution  is  obtained  by  an  elementary  one-dimensional 
calculation  exploiting  the  radial  symmetry.  This  calculation  uses  the  random  choice 
method,  together  with  a  time-splitting  step  that  takes  account  of  the  source  term  mtro- 
duced  by  the  change  in  coordinates  [22].  In  the  front  tracking  method  it  is  necessary  to 
track  the  contact  wave  as  well  as  the  shock,  since  otherwise  the  Lax-Wendroff  interior 
scheme  is  too  unstable  on  reasonable  grids. 

In  Figs.  4a-d  are  shown  the  results  of  a  calculation  on  a  40  by  40  grid  in  which  the 
pressure  inside  the  initial  circle  is  100  times  higher  than  outside.  Fig.  4a  shows  the  posi- 
tions of  the  contact  (the  inside  circle)  and  the  shock  (the  outside  circle)  together  with  the 
velocity  vectors  at  a  time  midway  in  the  run.  The  corresponding  density  and  pressure 
profiles  as  a  function  of  radius  are  shown  in  Figs.  4b,c.  In  these  profiles  the  solid  curve  is 
the  result  from  the  one-dimensional  calculation,  while  the  vertical  error  bars  extend  from 
the  minimum  to  the  maximum  values  in  the  two-dimensional  solution  at  fixed  radii,  indi- 
cating the  angular  dependence  of  the  solution.  In  Fig.  4d  the  contact  and  shock  positions 
are  plotted  vs.  time,  with  the  solid  curve  showing  the  comparison  answer,  and  the  dots 
showing  the  radii  of  the  two-dimensional  fronts.  The  deviation  of  these  fronts  from  cir- 
cles is  too  small  to  be  plotted  in  this  figure. 


-22- 

Similarly,  Rgs.  5a-d  show  the  results  for  a  calculation  in  which  the  pr  assure  on  the 
inside  is  100  times  smaller  than  on  the  outside.  Since  for  this  run  there  is  iin  outward- 
moving  rarefaction,  a  circular  reflecting  wall  (the  outer-most  circle  in  Rg.  5a)  was  intro- 
duced in  order  to  preserve  the  circular  symmetry  of  the  problem. 

6.2.  Snpersonic  Flow  Past  a  Wedge 

The  supersonic  flow  past  a  wedge  in  a  channel  was  tested  by  comparing  the  solution 
with  the  solution  of  the  steady-state  solution  obtained  by  the  method  of  characteristics 
[17].  In  this  test  a  bow  shock  generated  by  a  Mach  3  flow  over  a  30^  wedge  interacts 
with  a  Prandti-Meyer  expansion.  The  initial  data  for  the  two-dimensional  calculation  was 
a  shght  perturbation  of  the  steady  state  solution  obtained  Ujing  the  method  of  characteris- 
tics. The  flow  was  simulated  using  a  50  by  50  grid  for  200  time  steps,  so  that  an 
upstream  signal  had  time  to  move  across  the  computational  region. 

In  Rgs.  6a  and  6b  we  show  the  initial  and  final  shock  positions  together  with  the  iso- 
pycnic  (constant  density)  contours.  In  Figs.  7a  and  7b  we  show  the  initial  and  final  den- 
sity distributions  along  the  two  sides  of  the  shock.  Rgs.  8a,b  and  9a,b  show  the  emalo- 
gous  density  distributions  along  the  wedge  and  along  the  portion  of  the  exit  below  the 
shock,  respectively.  These  figures  indicate  that  the  front  tracking  scheme  accurately 
reproduces  the  steady  state  solution.  In  this  test  it  was  especially  important  to  apply  the 
boundary  conditions  along  the  wedge  correctiy. 

6.3.  The  Kelvin-Helmholtz  Instability 

The  classical  Kelvin-Helmholtz  instability  concerns  two  fluids  separated  by  an  inter- 
face across  which  there  is  a  discontinuity  in  the  tangential  velocity.  Such  a  flow  confi- 
guration is  imstable  against  a  sinusoidal  perturbation  of  the  interface.  In  the  regime 
where  the  ampUtude  is  small  relative  to  the  wavelength,  a  first  order  correction  to  the 
linearized  equations  provides  a  comparison  solution.  As  illustrated  in  Rg.  10,  we  look  for 
a  periodic  solution  of  the  equations  of  gas  dynamics  that  reduces,  in  the  Umit  where  the 
ampUtude  of  the  shp  line  is  small,  to  a  fiow  with  constant  density  and  pressure 
throughout,  and  constant  particle  velocity  vq  above  and  -vq  below.    By  linearizing  the 
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conservation  of  mass  and  Cnxxo's  eqii^itions  one  obtains  [19]  the  wave  equation  for  the 
velocity  potential.  The  solution  exhibit-s  an  exponentially  growing  sinusoidal  interface 
with  a  growth  rate  that  depends  on  the  jump  2  vq/cq  in  the  Mach  number  of  the  free 
stream  and  on  the  distance  to  the  boundaries.  The  corresponding  result  obtained  assum- 
ing incompressibility  of  the  fluid  is  obtained  in  the  limit  of  vanishing  Mach  number. 

The  results  of  our  numerical  experiments  are  shown  in  Figs.  ll,12a-c.  Fig.  11  is  a 
plot  of  the  amplitude  of  the  sinusoidal  slip  line  as  a  function  of  time  for  a  small  initial  per- 
turbation (amplitude-to-wavelength  ratio  of  0.005).  There  are  two  cases,  corresponding 
to  Mach  number  jump  1  and  0.1.  The  results  from  the  linear  analysis,  both  compressible 
and  incompressible,  are  superimposed.  As  this  plot  demonstrates,  the  growth  rate  is  accu- 
rately calculated  even  in  the  regime  where  methods  for  incompressible  fluids  are  invalid. 

Calculations  were  also  performed  in  the  large  amplitude  regime.  The  initial  flow 
field,  shown  in  Fig.  10,  was  obtained  using  the  perturbation  formulae  at  an  amplitude-to- 
wavelength  ratio  of  0.2  and  Mach  number  jump  0.4.  Rg.  12a-b  show  the  slip  line  and 
the  momentum  density  vectors  as  computed  with  a  20  by  20  grid  and  a  40  by  40  grid 
respectively.  Fig.  12c  shows  the  slip  line  for  the  same  run  as  computed  with  an  80  by  80 
grid.  The  slip  lines  for  aU  three  runs  are  superposed  in  Rg.  12d.  All  three  runs  agree  in 
the  general  form  of  the  slip  line.  The  finer  grids  show  evidence  of  vortex  formation  and 
roll-up  at  shorter  wavelengths.  We  stress,  though,  that  even  the  qualitative  shape  of  the 
slip  line  is  extremely  sensitive  to  the  initial  flow  conditions.  For  instance  if  the  initial  slip 
line  is  taken  to  be  purely  sinusoidal  in  the  vortidty  rather  than  in  the  horizontal  displace- 
ment [15],  the  slip  line  we  calculate  evolves  into  a  single  tight  spiral.  The  authors  know 
of  no  comparison  solution  for  this  problem,  which  involves  the  effects  of  the  compressibil- 
ity of  the  fluid  and  the  presence  of  boundaries,  but  the  results  are  in  qualitative  agree- 
ment with  incompressible  calculations. 

6.4.  Regular  Reflection  of  a  Shock  Wave 

Finally,  the  numerical  solution  for  non-steady  regular  reflection  of  a  shock  wave 
reflections  were  compared  with  experimental  results  [4].  In  the  test  problem  a  shock  with 
Mach  number  is  M^  =  2.05  is  incident  on  a  wedge  with  angle  ©vv  =  63.4°. 
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We  started  the  calculations  with  a  small  reflected  shcx±  enclosing  a  region  (called  the 
bubble)  one  quarter  of  a  mesh  interval  in  height  and  with  almost  arbitrary  data  inside. 
The  calculations  were  somewhat  unstable  in  the  neighborhood  of  the  comer  for  a  bubble 
smaller  than  two  mesh  intervals  in  height.  After  the  bubble  was  higher  than  four  mesh 
intervals  the  effects  of  the  instability  were  not  significant.  At  this  stage  of  the  run  wc 
found  that  the  solution  had  settled  down  to  its  self-similar  form  and  was  largely  indepen- 
dent of  the  initial  dats.  In  Fig.  13a  the  isopycnic  (constant  density)  contours  that  we 
obtained  numerically  arc  shown  at  two  stages  in  the  run.  In  Fig.  13b  the  wall  density  dis- 
tribution obtained  in  our  calculations  at  the  end  of  the  run  is  superimposed  on  the  experi- 
mental data. 

7.   Conclusions 

We  have  seen  that  front  tracking  gives  correct  results  for  a  series  of  test  problems  in 
compressible  fluid  dynamics.  Computations  on  even  coarse  grids  give  satisfactory  ar.swers 
for  the  problems  considered  here.  For  example,  in  the  regular  reOection  run,  convergence 
to  the  self-similar  solution  has  occurred  when  the  region  enclosed  by  the  reflected  shock  is 
about  five  mesh  blocks  high.  In  the  runs  with  circularly  symmetric  initial  data,  the  con- 
tact and  shock  are  initially  separated  by  0.4  mesh  blocks,  while  by  the  middle  of  the  rui5 
they  are  separated  by  3  mesh  blocks.  Successful  solution  of  test  problems  on  coarse  grids 
is  critically  important  if  complicated  problems  are  to  be  solved  on  computers  available 
within  the  near  future.  Also  the  ability  to  handle  bifurcations  in  the  front  topology  is 
important  for  complicated  problems.  This  capability  exists  at  present  for  the  oil  reservoir 
version  of  front  tracking  [9]  and  is  currentiy  being  implemented  for  gas  dynamics. 

Exactly  the  same  source  code  for  the  interface,  front,  hyperbolic,  driver,  and  utility 
libraries  have  been  applied  to  diverse  physical  applications.  The  front  tracking  code  is  a 
large  code,  but  its  modular  structure  has  minimized  the  problems  typically  associated  with 
the  development  and  use  of  such  codes.  The  multiple  applications  have  encouraged  a 
more  universal  approach  to  the  portions  of  the  code  for  which  generality  is  appropriate. 
Moreover,  each  application  has  benefited  from  progress  necessitated  by  special  difficulties 
in  other  applications.    See  [7,8,9,13]  for  validation  of  front  tracking  for  incompressible 
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flow  and  for  oil  reservoirs.    Thus  front  tracking  has  been  fruitful  in  providing  a  general 
framework  for  the  nun^rical  solution  of  hydrodynamic  problems. 
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Fig.  1.   The  local  coordinate  system  used  for  propag:itir3  iLc  front. 
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Fig.  2.  The  characteristic  curves  for  a  forward  shock. 


/    coM+ac't 


^  ^ 


Rg.  3.  The  characteristic  curves  for  a  contact  discontinuity. 
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Fig.  4a.  The  shock  wave  (outer  circk),  the  contact  discontinuity  (inner  circle),  and  the 
velocity  vectors  are  pictured  for  a  circularly  symmetric  computation  using  a  40  by  40  grid. 
The  initial  conditions  consisted  of  uniform  density,  zero  velocity,  and  a  circular  pressure 
discontinuity  at  radius  0.2  with  pjressurc  ratio,  inside  to  outside,  of  100.  The  time  it  takes 
a  sound  wave  in  the  region  inside  of  the  contact  to  travel  0.4  times  the  width  of  the  com- 
putational region  has  elapsed. 
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Fig.  4b.  A  plot  of  density  vs.  radius  corresponding  to  Fig.  4a  is  shown.  The"  solid  curve 
shows  the  results  obtained  in  a  one-dimensional  calculation  using  the  random  choice 
method.  The  vertical  lines  represent  the  range  of  density  values  in  the  two-dimensional 
solution  at  a  fixed  radius  as  the  angle  varies.  Thus  the  vertical  lines  show  the  angular 
dependence  in  the  solution. 
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Fig.  4c.  A  plot  of  pressure  vs.  radius  corresponding  to  Fig.  4a.  is  shown.  The  solid  curve 
and  the  vertical  lines  represent  the  one-  and  two-dimensional  results,  as  explained  in  the 
caption  to  Fig.  4b. 
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Fig.  Ad.  A  plot  of  the  radius  cf  the  contact  dkccntinuity  (Q  and  the  radius  of  the  shock 
wave  (S)  as  functions  of  tinie  h  presented.  The  soUd  curves  were  obtained  by  passive 
tracking  in  the  cnc-dinjcnsicnal  random  choice  solution.  The  dots  represent  radius  values 
in  the  two-dimensional  solution.  The  angular  dependence  of  the  radius,  i.e.  tl:^  deviation 
of  the  tracked  front  from  a  drelc,  is  too  small  to  plot. 


-33- 


Fig.  5a.  The  shock  wave  (inner  drclc),  the  contact  discontinuity  (outer  drcle),  and  the 
velocity  vectors  are  pictured  for  a  circularly  symmetric  computation  using  a  40  by  40  grid. 
The  initial  conditions  consisted  of  uniform  density,  zero  velocity,  and  a  circular  pressure 
discontinuity  st  radius  0.35  with  pressure  ratio,  inside  to  outside,  of  0.01.  The  time  it 
takes  a  sound  wave  in  the  region  outside  of  the  contact  to  travel  0.48  times  the  width  of 
the  computational  region  has  elapsed.  The  outermost  circle  is  a  reflecting  wall,  intro- 
duced to  mflintaiTi  circular  symmetry. 
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Fig.  5b.  A  plot  of  density  vs.  radius  corresponding  to  Fig.  5a  is  shown.  The  solid  curves 
and  the  vertical  lines  represent  the  one-  and  two-dimensional  results,  as  explained  in  the 
caption  to  Fig.  4b. 


-35- 


l.lf 


A-^x 


jij]Mi*^'^~^-!sxiii 


r«Ai 


us 


O.S 


Fig.  5c.  A  plot  of  pressure  vs.  radius  corresponding  to  Fig.  5a.  is  shown.  The  solid 
curves  and  the  vertical  lines  represent  the  one-  and  two-  dimensional  results,  as  explained 
in  the  caption  to  Fig.  4b. 
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Fig.  5d.  A  plot  of  the  radius  of  the  contact  discontimiity  (C)  and  the  radius  of  the  shock 
wave  (S)  as  functions  of  time  is  presented.  The  solid  curves  and  the  dots  represent  the 
one-  and  two-dimensional  solutions,  as  explained  in  the  caption  in  Fig.  4d. 
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Fig.  6a.  The  position  of  the  bow  shock  (S)  and  the  isopycnic  (constant  density)  contours 
are  shown  for  the  steady-state  flow  configuration  obtained  when  parallel  supersonic  flow 
with  Mach  number  3  impinges  from  the  left  on  a  30°  wedge  (W).  These  data  were  calcu- 
lated using  a  one-dimensional  random  choice  method  for  steady,  supersonic  flow. 
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Fig.  6b.  The  position  of  the  bow  shock  (S)  and  the  isopycnic  contours,  as  obtained  in  a 
time-dependent  calculation  on  a  50  by  50  grid  starting  from  the  steady-state  flow  condi- 
tions of  Fig.  6a,  arc  shown.  The  time  it  takes  a  sound  wave  in  the  region  ahead  of  the 
bow  shock  to  travel  upstream  1.2  times  the  length  of  the  wedge  has  elapsed. 
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Fig.  10.  The  flow  configuration  for  the  Kelvin-Helmholtz  instability  is  plotted.  There  is 
a  sinusoidal  slip  line  (solid  line)  separating  the  region  above  from  the  region  below.  In 
the  regime  whore  the  amplitude  of  the  slip  line  is  small  the  density  and  pressure  arc  al- 
most constant  throughout  the  flow,  while  the  fluid  velodtics  arc  approximately  horizontal, 
with  the  velocity  above  equal  and  opposite  to  to  the  velocity  below. 
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Fig.  11.  A  plot  of  the  amplitude  of  the  slip  line  vs.  time  is  shown  for  the  Kelvin- 
Ifelmholtz  instability  in  the  small  amplitude  regime.  For  the  curves  (L  j  -  L  3)  (which  are 
indistinguishable  in  this  graph)  the  jump  in  the  velocity  has  Mach  number  0.1,  while  for 
the  curves  (H  j  —  H  3)  the  jump  in  the  velocity  has  Mach  number  1.  In  each  case  curve  1 
was  obtained  numerically  using  a  20  by  20  grid,  whereas  curves  2  and  3  were  obtained  for 
compressible  and  incompressible  flow,  respectively,  using  perturbation  analysis  in  the  am- 
plitude of  the  slip  line.  At  the  higlter  Mach  number  it  is  important  to  model  the  signifi- 
cant compressibility  effects. 
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Fig.  12a.  The  slip  line  and  the  momentum  densities  are  shofwn  for  the  Kelvin-Hclmholtz 
instability  with  a  large  amplitude.  The  calculation  was  made  on  a  20  by  20  grid.  The  ini- 
tial data,  shown  in  Fig.  10,  was  obtained  by  the  perturbation  analysis  for  a  velocity  jump 
with  Mach  number  0.4.  The  time  it  takes  a  sound  wave  to  travel  0.72  wavelengths  has 
elapsed.  In  the  large  amplitude  regime  the  slip  line  no  longer  remains  sinusoidal,  but 
rather  rolls  up  as  vortices  are  formed. 
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Fig.  12b.   The  slip  line  and  the  momentum  densities  are  plotted  for  tlie  problem  described 
in  Fig.  12a,  as  calculated  on  a  40  by  40  grid. 
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Fig.  12c.   The  slip  line  position  is  plotted  for  the  problem  described  in  Fig.  12a,  as  calcu- 
lated on  a  80  by  80  grid. 
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Rg.  12d.  The  slip  line  positions,  as  computed  on  20  by  20,  40  by  40,  and  80  by  80  grids, 
are  superposed  for  the  problem  described  in  Fig.  12a.  Fine  structure  in  the  vortices  and 
in  the  shape  of  the  slip  line  becomes  evident  under  refinement  of  the  computational  grid. 
The  authors  know  of  no  comparison  solution  for  this  problem. 
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Fig.  13a.  The  isopycnic  (constant  density)  contours  are  shown  for  a  regular  reflection 
problem  in  which  an  incident  shock  (I),  moving  downward  and  to  the  right  into  ambient 
gas,  impinges  on  a  reflecting  wall  (W)  and  reflects  at  point  (r)  to  form  a  reflected  shock 
(R).  In  this  computation  the  Mach  number  of  the  incident  shock  is  2.05,  and  the  angle  of 
incidence  of  the  shock  on  the  wall  is  26.6°.  The  calculations  were  performed  on  a  60  by 
40  grid. 
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Fig.  Db.  The  density  distribution  along  the  wall  for  the  problem  described  in  Fig.  Da  is 
plotted  vs.  distance  along  the  wzdl.  The  solid  curve  was  obtained  by  nunaerical  calcdla- 
tion,  while  the  dots  show  the  experimental  results  for  air  obtained  by  Deschambault  and 
Glass. 
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